Phytochemistry reflects different evolutionary history in traditional classes versus specialized structural motifs

Foundational hypotheses addressing plant–insect codiversification and plant defense theory typically assume a macroevolutionary pattern whereby closely related plants have similar chemical profiles. However, numerous studies have documented variation in the degree of phytochemical trait lability, raising the possibility that phytochemical evolution is more nuanced than initially assumed. We utilize proton nuclear magnetic resonance (1H NMR) data, chemical classification, and double digest restriction-site associated DNA sequencing (ddRADseq) to resolve evolutionary relationships and characterize the evolution of secondary chemistry in the Neotropical plant clade Radula (Piper; Piperaceae). Sequencing data substantially improved phylogenetic resolution relative to past studies, and spectroscopic characterization revealed the presence of 35 metabolite classes. Metabolite classes displayed phylogenetic signal, whereas the crude 1H NMR spectra featured little evidence of phylogenetic signal in multivariate tests of chemical resonances. Evolutionary correlations were detected in two pairs of compound classes (flavonoids with chalcones; p-alkenyl phenols with kavalactones), where the gain or loss of a class was dependent on the other’s state. Overall, the evolution of secondary chemistry in Radula is characterized by strong phylogenetic signal of traditional compound classes and weak phylogenetic signal of specialized chemical motifs, consistent with both classic evolutionary hypotheses and recent examinations of phytochemical evolution in young lineages.

www.nature.com/scientificreports/ high levels of divergence among Piper clades. For Bayesian phylogenetic inference, we mitigated the influence of missing data by removing loci absent in > 30% of samples. The final dataset for phylogenetic analysis consisted of 641 ddRADseq loci (~ 86 bp in length each) that housed 9113 genetic variants (51% parsimony informative). Aligned loci were concatenated into a nexus alignment with missing data at 18.9% of sites. Bayesian phylogenetic analysis of ddRADseq data resolved eight major Neotropical Piper clades with high posterior support (Fig. 1). While past phylogenetic studies supported the monophyly of seven of these eight clades (Macrostachys, Radula, Peltobryon, Pothomorphe, Hemipodium, Isophyllon, and Schilleria) 34,43 , our analysis resolved an additional clade, Churumayu. Notably, Isophyllon and Churumayu were highly supported, monophyletic clades and not nested within Radula, as was inferred in previous analyses 43 . Contrary to previous phylogenetic hypotheses of Piper 34,43 , our analysis might suggest Churumayu is the most basal clade, but we caution that this node had very low posterior support (51%). Intrageneric relationships below the clade level were highly resolved, with nearly all nodes exhibiting greater than 95% posterior support, including within the diverse Radula clade (Fig. 1). Our phylogenetic hypothesis for Radula indicates three species (P. hispidum, P. colonense, P. lucigaudens) may be paraphyletic.
Phytochemical diversity in Piper. All but four individuals included in the inferred Piper tree were successfully chemically extracted and profiled. Nearly all common compound classes that have been previously reported in Piper 46 were observed from our compound characterization analysis (see Table S2). This analysis Red circles indicate support values between 75 and 85%. Three nodes with less than 75% posterior support were given numerical support values. Blue bars at each node denote the 95% highest posterior density interval on relative node ages. The photos to the right of the tree showcase a sample of Piper diversity, including a few of the species which were included in this study: (a) Piper hillianum (Macrostachys), (b) P. acutifolium (Peltobryon), (c) P. umbellatum (Pothomorphe), (d) P. pseudofuligineum (Radula), (e) P. concepcionis (Radula), (f) P. disparipes (Radula), (g) P. friedrichsthalii (Radula), (h) P. dilatatum (Radula), (i) P. bredemeyeri (Radula), (j) P. immutatum (Radula), (k) P. erubescentispicum (Radula), and (l) the widespread and often weedy P. aduncum (Radula www.nature.com/scientificreports/ revealed the presence of broad metabolite classes that are ubiquitous across plant families (e.g., lignans, flavonoids/chalcones, etc.) as well as classes that are specifically common in Piper (e.g., amides) (Fig. 2, Table S2). Specific compound characterization revealed genus specific compounds and compound classes (piplartine, cenocladamide, crassinervic acid, kava lactones), as well as metabolites that are more rarely reported in plants (putrescine diamides, nerolidyl catechol, alkenyl phenols, anuramide peptides) (Fig. 2, Table S2). Alternative methods, such as sampling across a species' ontogeny, sampling reproductive parts or roots, and storing freshly collected tissue in methanol rather than air drying would add to a more comprehensive picture of variation in phytochemical diversity across and within species, but our sampling was standardized to allow for initial comparisons across species, some of which were collected in remote regions.
Metabolite phylogenetic signal and evolutionary associations. We recovered 35 metabolite classes, of which only eight were sufficiently present across our taxa to afford tests of phylogenetic signal and correlated evolution. For all eight metabolite classes, estimates of D did not deviate from a null distribution generated under a scenario of Brownian motion (Table 1), consistent with phylogenetic signal. Two of the eight traits, phenolic glycosides and lignans, exhibited strong phylogenetic signal (D < 0), while the remaining six traits exhibited weak phylogenetic signal (0 < D < 1). Further, all but one of the metabolite classes had observed values of D that differed from a null distribution generated under a phylogenetic randomness scenario ( www.nature.com/scientificreports/ and (2) p-alkenyl phenols and kavalactones/butenolides. For the first pair of traits, a model of contingency in which changes in chalcones depend on the state of flavonoids provided the best fit to the data ( Table 2). In this model, when flavonoids are present, chalcone gains are 1.4 times more probable than chalcone losses; however, when flavonoids are absent, chalcone losses are much more probable than chalcone gains (Fig. 3). The alternative contingency model for this pair of traits (i.e., changes in flavonoids depend on the state of chalcones) was also a good fit to the data ( Table 2). According to this model, when chalcones are present, flavonoid gains are approximately nine times more probable than flavonoid losses. Alternatively, when chalcones are absent, flavonoid losses are approximately five times more probable than flavonoid gains (Fig. 3). For the second pair of traits, p-alkenyl phenols and kavalactones/butenolides, the best fit model was one of interdependent evolution in which changes in p-alkenyl phenol depend on the state of kavalactones/butenolides, and vice versa ( Table 2). When kavalactones/butenolides are present, p-alkenyl phenol transitions are more probable than when they are absent, with the loss of p-alkenyl phenols being much more probable than the gain of p-alkenyl phenols under both scenarios. Alternatively, when p-alkenyl phenols are present, the loss of kavalactones/butenolides is extremely probable relative to the gain of kavalactones/butenolides, which is rarely observed. When p-alkenyl phenols are absent, kavalactones/butenolides are rarely gained or lost (Fig. 3).
Phylogenetic signal in high-dimensional metabolomic data. While the eight metabolite classes uniformly exhibited at least moderate levels of phylogenetic signal, evidence for phylogenetic signal in multivariate analyses of the crude 1 H NMR data was largely absent. PCo axes 1 & 2 and 3 & 4 explained 32.8% and 16.0% of variance in the 1 H NMR data, respectively, but showed little clustering by clade (Fig. 4a). Permutational multivariate analyses of variance were not significant for combinations of either PCo 1 & 2 (P = 0.407) nor 3 & 4 (P = 0.142), suggesting that different clades do not form distinct clusters in chemospace based on their 1 H NMR spectra.
According to the MRM models, phylogenetic distance significantly predicts phytochemical distance within Radula (β = 4.503, P = 0.013) but not across all clades (β = 1.775, P = 0.146) (Fig. 4b). It is important to note that the Table 1. Estimates of phylogenetic signal (D) 57 Table 2. Correlated evolution was detected in two pairs of metabolite classes with Pagel's method 76 : (1) chalcones and flavonoids; and (2) kavalactones/butenolides and p-alkenyl phenols. A model comparison framework was employed to evaluate four potential models of trait evolution using AIC: interdependent evolution (transition rate in one trait depends on state at another, and vice versa); contingent change (transition rate in one trait depends on state at another, but not the converse); and independent evolution.  www.nature.com/scientificreports/ proportion of variance explained by the significant MRM model is low (R 2 = 0.039), suggesting that the majority of variation in NMR data cannot be explained by phylogenetic distance. Analyses with the generalized K statistic (K mult ) indicated lower levels of phylogenetic signal in the metabolomic data than expected under a Brownian motion model of evolution for Piper generally (K mult = 0.1606, P = 0.001) and for Radula specifically (K mult = 0.1803, P = 0.001). Still, the observed K mult was higher than all K mult values obtained with permutations of the 1 H NMR dataset (Fig. S1). Additionally, few K mult tests of the permuted data yielded significant P-values (4.4% of permutations), indicating that the estimate we observed, though subtle and lower than Brownian motion expectations, was real and not a statistical artifact of zero-inflation in the data.

Discussion
Piper is a hyper-diverse lineage in which phytochemical diversity has influenced evolutionary and ecological processes and shaped complex tropical communities 15,39 . However, limitations in both the degree of phylogenetic resolution and the understanding of phytochemical diversity in this group have precluded analyses of phylogenetic signal and correlated evolution of phytochemistry. Phylogenies inferred here with ddRADseq data substantially improved resolution and support compared to past studies of Piper, which were limited by interspecific variation in small numbers of Sanger-sequenced loci 34,42,43 . Although the data set did not include members from all previously recognized groups, analyses resolved eight monophyletic Neotropical Piper clades, seven of which have been inferred in previous analyses of the genus based on chloroplast psbJ-petA and nrITS 34,43 . Two of the eight clades, Churumayu and Isophyllon, had been previously nested within Radula 43 ; however, our results suggest that they are independent monophyletic lineages (Fig. 1). Despite low support for several deep divergences, the phylogeny inferred here had strong resolution and support for recent relationships, including within Radula (Fig. 1), consistent with other recent reduced representation sequencing studies that have generated high quality phylogenies at shallow time scales 28,31,32 . However, a potential limitation of such sequencing designs may include the recovery of fewer loci shared by more distantly related samples due to allelic dropout 47 . It is possible that allelic dropout, potentially exacerbated by strict filtering based on missing data, led to weak support values for deep splits in the phylogeny, many of which occurred early in the history of the Neotropical  www.nature.com/scientificreports/ Piper lineage 34 . Nonetheless, the resulting subset of data (641 loci; 9113 SNPs) was sufficient for inferring a largely resolved phylogeny, highlighting the potential promise of reduced representation sequencing for resolving evolutionary histories even in groups spanning moderately deep divergence. Although our sampling was limited to 44 of 450 estimated species within Radula, the extent of sampling is a substantial improvement over past phylogenetic analyses for the group 42,43 .
Comparative studies have taken diverse approaches to analyzing metabolomic data, each providing a unique perspective on the evolution of specialized metabolites 10,24 . Here, we first characterized the presence/absence of 35 metabolite classes commonly used to categorize plant secondary compounds that are hierarchically nested into three levels of structural resolution. Specific categories at the lowest level of the hierarchy, representing specialized structural motifs or specific molecules, were rare across species and precluded tests of phylogenetic signal and correlated evolution at our level of taxonomic sampling (Fig. 2). Despite not being able to test for phylogenetic signal, clustering is evident for more specific categories, such as crassinervic acid and prenylated flavonoids, which are only present in small subclades but include particularly effective defenses 36,46 . Alternatively, broader metabolite classes at intermediate and high positions in the hierarchy that are directly tied to fundamental secondary metabolite biosynthetic pathways were more abundant across species and exhibited moderately high levels of phylogenetic signal across Radula (Table 1, Fig. 2). This pattern may be expected if initial biosynthetic steps are conserved over longer evolutionary scales, permitting the abundance of broad chemical classes, yet later stage modifications of these core structures are more evolutionarily labile, causing structural similarity to be low even among related species. Flavonoids are a good example of this pattern, with pathways that form the flavonoid scaffold being very conserved, as they are catalyzed by modified enzymes from ubiquitous metabolic pathways, but then subsequent biosynthetic steps (e.g., those catalyzed by p450 enzymes) modify these scaffolds 48 , yielding unique molecules towards the tips of evolutionary trees (Fig. 3E). For example, late-stage modification of common flavonoid scaffolds can result in the production of non-aromatic protoflavonoids. These compounds rarely occur across the plant kingdom and have only recently been found in one species of Piper 49 , but this type of subtle structural modification that leaves most of the flavonoid scaffold intact dramatically enhances the cytotoxic properties compared to that of the parent flavonoid 50,51 .
One key prediction from the escape and radiate hypothesis is that adaptive defensive traits should be phylogenetically conserved within the lineage they evolved, but this prediction has mostly been evaluated with broad classes of secondary metabolites at high taxonomic scales 6,13,48 rather than specific compounds in recent diversifications 7,10,16 . A growing number of studies conducted at shallow evolutionary scales suggests low phylogenetic signal in many chemical traits 14,15,18 . While evidence for low phylogenetic signal is often attributed to high evolutionary rates (i.e., evolutionary lability), simulations under various evolutionary processes and conditions indicate that the relationship between phylogenetic signal and rate of trait evolution is not necessarily straightforward, and evidence for low phylogenetic signal is not an indication of any single evolutionary process 52 . Nonetheless, understanding how phylogenetic signal responds to variation in phylogenetic scale is informative in a comparative sense, especially among different traits or classes of traits generated with different levels of analytical resolution. Phylogenetic signal is also a useful starting point for developing insights into the drivers of herbivorous insect radiations, as codiversification in many of these lineages is structured in part by chemical defense and biotic interactions 40,53 . Our results are generally consistent with the predictions of moderately strong signal for broad classes of compounds, as well as the lack of signal for specific structures captured by 1 H NMR data.
The 1 H NMR data address a different set of hypotheses than data from categorization of individual molecules-peaks represent resonances associated with particular molecular structures rather than individual compounds, and the chemical shift (frequency), shape, and abundance of these resonances are extremely sensitive to subtle structural changes. 1 H NMR spectroscopy easily detects a great range and subtle differences in compositional and structural complexity, including increasing size, asymmetry and oxidation states, that might be predicted to evolve in response to divergent selection across plant populations responding to different suites of enemies 22 . Low levels of phylogenetic signal in the 1 H NMR data is also likely due to the fact that many molecular features of small defensive molecules have potentially evolved in a convergent manner across Piper, such as the kavalactones, p-alkenyl phenols, piplartine, oxidized prenylated benzoic acids, chromenes, anuramide peptides, and phenethyl amides.
There are numerous limitations that could affect estimates of phylogenetic signal in comparative studies 54 that are relevant to the analyses presented here. First, incomplete taxon sampling likely influenced our results to some degree, but sampling was conducted randomly, and the probability that a particular species was sampled was unlikely related to any aspect of its chemical phenotype 55 . Low sampling proportion in clades other than Radula may have reduced our power to detect phylogenetic signal across all our sampled clades 55 (Fig. 4b). However, despite only sampling approximately 10% of the Radula clade of Piper, our sample size should provide sufficient power to infer phylogenetic signal in this clade if present 56,57 (Fig. 4b). Second, while topological errors and small sample size may have reduced our power to detect phylogenetic signal at deeper time scales 58 , more comprehensive genomic sampling produced enhanced phylogenetic resolution of the Radula clade, where we focused the majority of phylogenetic comparative methods. In addition, we were unable to quantify the measurement error associated with the chemical traits within species, which can decrease the statistical power for detecting phylogenetic signal 56,59,60 . It is also possible that environmental effects on our chemical traits could bias estimates of phylogenetic signal and correlations 59 .
The causes of correlated evolution, including linkage, epistasis, and selection, are difficult to detect without careful approaches in quantitative genetics and population genomics. Nevertheless, one advantage of examining the presence/absence of multiple classes of defensive compounds in a phylogenetic context is that it is possible to test for expected patterns of correlated evolution due to shared metabolic pathways (e.g., flavonoids and cardenolides 7 ) or due to adaptive advantages of specific mixtures. Recent studies detecting evolutionary www.nature.com/scientificreports/ associations among chemical traits 17,18 have posited that the branching structure of metabolic pathways could potentially drive this pattern. If metabolite classes share a common precursor, one might expect evolutionary tradeoffs and negative covariation. Alternatively, if metabolite classes lie along the same metabolic pathway, an increase in one class may be concomitant with increases in another (or vice versa), causing positive covariation among the classes. There are also numerous empirical examples supporting the hypotheses that correlations may be driven by functional redundancy 61 or selection for synergistic effects on herbivores 20 rather than the structural constraints of metabolism. Suites of covarying defensive traits, or defense syndromes, have been detected in several plant genera 9,53 and plant communities 62 , and have been predominantly used to describe covariation among mechanical and chemical defenses. It is interesting to note the correlated evolution of the flavones/chalcones and the p-alkenyl phenols/kavalactones could be due to metabolic constraints, as well as possible adaptations via synergistic (e.g., kavalactones in P. methysticum) or other mixture-associated defensive attributes 22 . Flavonoids and chalcones are directly linked biosynthetically, such that the inherent reactivity of the chalcone moiety permits the enzymatic processes that result in cyclization to the flavonoid scaffold (Fig. 3e). This strong biosynthetic tie yields a clear prediction that the presence of one would depend on the other, and indeed our structural analysis found many cases where both metabolite classes co-occurred in the same sample. Revealing the relationship between the kavalactones and p-alkenyl phenols is more tenuous because both classes are less prevalent across our samples. Kavalactones and p-alkenyl phenols are dramatically different compounds that diverge at a much earlier branch point from a common cinnamic/coumaric acid precursor. Whereas one polyacetate chain extension pathway leads to the long-chain lipophilic substituent, characteristic of the p-alkenyl phenols, the other chain extension pathway conserves oxidation states through the chain extension process to produce the lactones (kavalactones or butenolides) through cyclization reactions (Fig. 3e). The overall outcome is different than the chalcone-flavonoid relationship; in this case, two dramatically different compounds are produced by divergence from a common early-stage biosynthetic precursor in contrast to the immediate biosynthetic precursor relationship between chalcones and flavonoids. Broader sampling across Piper and Radula will be necessary to confirm this unexpected relationship between kavalactones and p-alkenyl phenols.

Conclusion
Here we sought to advance understanding of phylogenetic relationships within Piper while simultaneously investigating the mode and manner of phytochemical evolution in this group. In addition to generating a well-resolved phylogeny, our results support theoretical expectations that broad classes of compounds display higher degrees of phylogenetic signal than molecular features revealed by 1 H NMR data. In addition, trait associations observed in Radula can be used to pose functional hypotheses about genetic constraints or biases on phytochemical evolution and how these factors structure plant-animal interactions. Such investigations are one of the emerging frontiers in terrestrial ecology, and we hope that our study provides one example of how collaborative and multi-disciplinary research can progress in this area. T. in the field, and confirmed with vouchers in the herbarium using regional keys, where available, comparison with type specimens, and experience with the genus. For chemical profiling and DNA sequencing, we collected the youngest, fully expanded leaves and dried them immediately with silica gel. While drying on silica gel may not inhibit enzymatic activity and could limit our analyses to relatively stable molecules, this is not an issue for the phylogenetic analyses described below. Collections were only made from mature individuals in the field. Vouchers were pressed, dried, and deposited in one or more herbaria for future reference and species verification (Table S1). To investigate the evolution of phytochemistry at a relatively shallow evolutionary scale, we conducted the majority of our sampling within Radula 34 .

Methods
Phylogenetic analyses. Genome-wide polymorphism data was generated for 71 individuals for phylogenetic analyses. Either the same accession sampled for chemical analysis, or an individual from the same population as the one sampled, were sequenced with a genotyping-by-sequencing approach 63 that is analogous to ddRADseq 64 . Briefly, genomic DNA was digested with two restriction enzymes, EcoRI and MseI. Sample-specific barcoded oligos containing Illumina adaptors were annealed to the EcoRI cut sites, and oligos containing the alternative Illumina adaptor were annealed to the MseI cut sites. Fragments were PCR amplified and pooled for sequencing. The library was size-selected for fragments between 350 and 450 base pairs (bp) with the Pippin Prep System (Sage Sciences, Beverly, MA), and sequenced on two lanes of an Illumina HiSeq 2500 at the University of Texas Genome Sequencing and Analysis Facility (Austin, TX). Single-end, 100 bp, raw sequence data were filtered for contaminants (E. coli, PhiX, Illumina adaptors or primers) and low quality reads using www.nature.com/scientificreports/ bowtie2_db 65 and a pipeline of bash and perl scripts (https:// github. com/ ncgr/ tapio ca). We used custom perl scripts to demultiplex our reads by individual and trim barcodes and restriction site-associated bases. Assembly and initial filtering was conducted with ipyRAD v.0.7.30 66 . ipyRAD was specifically designed to assemble ddRADseq data for phylogenetic applications, permits customization of clustering and filtering, and allows for indel variation among samples 66 . Because a suitable Piper genome was not available at the time of analysis, we generated a de novo consensus reference of sampled genomic regions with ipyRAD. Briefly, nucleotide sites with phred quality scores lower than 33 were treated as missing data. Sequences were clustered within individuals according to an 85% similarity threshold with vsearch 67 and aligned with muscle 68 to produce stacks of highly similar ddRADseq reads (hereafter, ddRADseq loci). The sequencing error rate and heterozygosity were jointly estimated for all ddRADseq loci with a depth > 6, and these parameters informed statistical base calls according to a binomial model. Consensus sequences for each individual in the assembly were clustered once more, this time across individuals, and discarded if possessing > 8 indels (max_Indels_ locus), > 50% heterozygous sites (max_shared_Hs_locus), or > 20% variable sites (max_SNPs_locus). To reduce the amount of missing data in our alignment matrix, ddRADseq loci were retained if they were present in at least 50 of 71 samples. The nexus file of concatenated consensus sequences for each individual, including invariant sites, was used as input for the Bayesian phylogenetic methods described below. Individual FASTQ files, nexus alignment, and complete information on additional parameter settings for this analysis are archived at Dryad (https:// doi. org/ 10. 5061/ dryad. j6q57 3nc7).
To resolve patterns of diversification and to provide a foundation for investigating variation in patterns of phytochemical evolution, we estimated a rooted, calibrated tree according to a relaxed clock model in RevBayes v.1.0.12 69 , which provides the ability to specify custom phylogenetic models for improved flexibility compared with other Bayesian approaches. The prior distribution on node ages was defined by a birth-death process in which the hyper priors on speciation and extinction rates were exponentially distributed with λ = 10. We relaxed the assumption of a global molecular clock by allowing each branch-rate variable to be drawn from a lognormal distribution. After comparing the relative fits of JC, HKY, GTR, and GTR + Gamma nucleotide substitution models with Bayes factors, we modeled DNA sequence evolution according to the best-fit HKY model. Eight independent MCMC chains were run for 100,000 generations with a burn-in of 1,000 generations and sampled every 10 generations. Chains were visually assessed for convergence with Tracer v.1.7.1 70 and numerically assessed with effective sample sizes (ESS), the Gelman − Rubin convergence diagnostic 71 , and by comparing the posterior probabilities of clades sampled between MCMC chains. The maximum clade credibility (MCC) tree provided the ultrametric fixed tree topology and relative node ages for phylogenetic comparative methods described below.
Chemical profiling. Crude proton nuclear magnetic resonance ( 1 H NMR) spectroscopy was chosen for chemotype mapping due to its ability to characterize subtle structural variation across a wide range of compound classes in a single, reproducible, non-destructive analysis 39 . Briefly, after leaf samples were ground to fine powder, approximately 100.0-2000.0 mg of leaf material were ground and transferred to a glass screw cap test tube with 10 ml of methanol, sonicated for 10 min, and filtered. This step was repeated and both filtrates were combined in a pre-weighed 20 ml scintillation vial. The solvent was removed in vacuo and dissolved in 0.6 ml methanol-d 4 for 1 H NMR analysis. Crude 1 H NMR solutions were standardized to 13.1 ± 3.8 mg/mL when possible and analyzed on a Varian 400 MHz solution state NMR spectrometer with autosampler. Data were processed using MestReNova software (Mestrelab Research, Santiago de Compostela, Spain). Spectra from the crude extracts were aligned with the solvent peak (CD 3 , δ = 3.31 ppm), baseline corrected, and phase corrected. Solvent and water peaks were removed and the binned spectra were normalized to a total area of 100. This data set is referred to as "crude 1 H NMR".
In addition to crude 1 H NMR spectral chemotyping, we further annotated samples based upon the presence or absence of compound classes. To further gain structural resolution across the crude extracts that were sampled, aliquots of the 1 H NMR extracts were diluted and subjected to GC-MS and LC-MS analysis (see Supplementary Information for additional details). Crude extracts were classified using chemotaxonomic classifications outlined in Parmar's comprehensive review of Piper phytochemistry 35 , and our rationale for assigning chemical classes is outlined for each species in Table S2. Briefly, phenolic compounds were identified from high-resolution matches to the METLIN mass spectrometry database 72 . Database hits were then confirmed by agreement of crude 1 H NMR chemical shifts with literature values for phenolics known to be found in Piper, but not always Radula species. Many compounds identified by LC-MS as flavonoids and chalcones had multiple possible METLIN matches, which confounded NMR confirmation. In these cases, we were still able to differentiate flavonoids from chalcones by characteristic UV spectra (l max ~ 350 nm). Phenylpropanoids and p-alkenyl phenols were identified based on characteristic GC-MS fragmentation for these compound classes known to be found in Piper. Piper amides were characterized in a similar fashion, starting from high-resolution mass spectrometric matches and confirming with known 1 H NMR data from the literature. In some cases, crude 2D-NMR analysis (COSY, HSQC) was used to confirm structural classifications. COrrelated SpectroscopY (COSY) was used to identify 1 H NMR that were contained within the same molecule, while Heteronuclear Single Quantum Coherence (HSQC) spectroscopy was used to identify the carbon ( 13 C) resonances associated with certain proton ( 1 H) signals to verify the presence of specific functional groups 73 . Only the most abundant and spectroscopically apparent compounds were classified due to the low sensitivity of NMR. 35  www.nature.com/scientificreports/ was not characterized for the sesquiterpene class. This hierarchical set of 35 traits is referred to as "metabolite classes" (Fig. 2). Additional details on chemical profiling can be found in the Supplementary Information.

Phylogenetic signal and evolution of metabolite classes.
To assess whether metabolite classes were phylogenetically conserved across Radula, we quantified phylogenetic signal in these binary traits using the D statistic 57 . The D statistic calculates the sum of sister-clade differences, Σd obs for an observed tree and binary trait, and scales this value with the distributions of sums expected under two disparate evolutionary models, random and Brownian motion (Σd r and Σd b , respectively), using the following equation: Thus, D is expected to equal 1 when the observed binary trait is distributed randomly, lacking phylogenetic signal, and is expected to equal 0 when it exhibits phylogenetic signal as expected under Brownian motion. As tests of phylogenetic signal with the D statistic are most accurate when the ratio of presences and absences is closer to 1:1 57 , we tested for phylogenetic signal in eight of the 35 metabolite classes (outlined in white in Fig. 2) which were present in a sufficient proportion of taxa. We used the phylo.d function in the caper package 74 in R v.4.0.0 75 to calculate the observed D for a subset of binary traits that were sufficiently present across the phylogeny. This value was compared to a distribution of D values simulated under models of phylogenetic randomness (D = 1) and pure Brownian motion (D = 0) to determine whether the observed D differed from either zero or one.
To detect evolutionary associations among pairs of metabolite classes within Radula, we used Pagel's method 76 that models evolutionary changes in two binary traits, X and Y, as continuous-time Markov processes in which the probabilities of state transition at one trait may depend on the state at the other trait. We tested all pairwise associations among the eight metabolite classes that were represented by a sufficient number of Radula taxa to provide accurate tests of evolutionary associations (N = 28). Significant tests of correlated evolution were followed by tests of contingency, in which changes at X depend on the state of Y, or vice versa. Model fits, comparisons, and plots were performed with the fitPagel function in the phytools package 77 in R.
Multivariate analyses of phylogenetic signal with crude 1 H NMR spectra. While the analyses above based on broad classifications of structurally determined metabolites provide a coarse view of phytochemical evolution, these classifications are anchored to the foundations of plant secondary metabolite biosynthesis. Using 1 H NMR spectra as a raw chemotype should allow a more detailed multivariate perspective on phytochemical diversity. Studies on other plant taxa have typically detected some signal and evolutionary correlations for broad classes of compounds but not necessarily for specific compounds or biologically active moieties, both of which can be inferred from 1 H NMR data. Multivariate approaches to phylogenetic comparative methods have provided insight into covarying suites of related traits, while simultaneously increasing the statistical power to detect phylogenetic signal 78 and differences in trait means among taxa 79 . Indeed, these multivariate approaches might be particularly useful when exploring the evolution of complex phenotypes, like the plant metabolome, which exhibit trait covariances due to metabolomic or functional associations 20 . Here we utilize three multivariate methods to detect patterns of phylogenetic signal for 263 resonances found in the crude 1 H NMR data representing all 35 metabolite classes: (1) principal coordinate analyses (PCoA); (2) multiple regression on distance matrices (MRM); and (3) multivariate estimation of phylogenetic signal.
To visualize patterns of chemotypic variation across all sampled species from all clades, we first analyzed the 1 H NMR data with PCoA. First we calculated the Manhattan distances between all pairwise species with the dist function in R, and then conducted PCoA on the distance matrix using the pcoa function in R. If the major axes of metabolomic variation are phylogenetically conserved, the plotted species scores should be clustered by clade in a rotated principal coordinate (PCo) space. Alternatively, if metabolomic variation is randomly distributed across the phylogeny, there should be little to no clustering by clade 80 . The degree to which plant clade predicted chemical similarity was assessed using permutational multivariate analysis of variance (permanova) 81 in the vegan package 82 in R based on Euclidean distances of the first four PCo axes.
Mantel tests have been frequently used to assess the degree of phylogenetic signal in multivariate data 10,83,84 by estimating the relationship between phylogenetic and phenotypic distances. Simulations under scenarios of measurement error have found instances where Mantel tests outperform traditional univariate methods in detecting phylogenetic signal, especially as the number of traits increases 60 . Because we were unable to account for measurement error in our study, we utilized Multiple Regression on distance Matrices (MRM) 85 to examine the relationship between metabolomic and phylogenetic distance at two evolutionary scales (within Radula and across all clades). Euclidean distances were calculated from the crude 1 H NMR spectra using the dist function in R, and phylogenetic distances for Radula only and all clades were calculated using the cophenetic.phylo function in the ape package 86 in R. MRM analyses were implemented using the MRM function with 1000 permutations in the ecodist package 87 in R.
Since Blomberg's K 56 statistic exhibits higher statistical power to detect phylogenetic signal relative to Mantel tests 88 , we quantified phylogenetic signal of the crude 1 H NMR at both evolutionary scales using a multivariate generalization of the K statistic (K mult ) 89 with the physignal function in the geomorph package 90 in R. Similar to the aforementioned D statistic, the K statistic compares the observed variation to that expected under Brownian motion, but the K statistic does not scale this comparison by the variation exhibited under a completely random evolutionary model 56,89 . Values of K greater than 1 indicate phylogenetic signal greater than expected under Brownian motion, whereas values between 0 and 1 indicate less signal than expected under Brownian motion. Significance for the generalized K statistic was assessed by permuting the 1 H NMR peak data among the tips of the phylogeny for 999 iterations. To determine whether the zero-inflated nature of the 1 H NMR data influenced the www.nature.com/scientificreports/ detection of phylogenetic signal, we permuted our 1 H NMR data set over 1000 iterations by randomly indexing our original 1 H NMR data matrix. This permutation method preserves the original proportion of zeros in the matrix while obfuscating any observed phylogenetic signal. The generalized K statistic test was calculated for each permutation, and our observed generalized K statistic was compared to the null distribution of permuted values.